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The discussion began with Dr. Bill Gropp introducing the concept of visual artifacts in 
numerical solutions. He presented examples of errors that appeared to be significant to the 
human eye, but that were well below the error criteria for the problem, and did not impact 
the quality of the numerical solution. 

The discussion then focused on defining a model problem where visual artifacts cduld be 
examined explicitly. 

Model Problem ^9^-417 

A. Louise Perkins _ ? 

Massachusetts Institute of Technology 

1 Introduction 

At the ICASE Workshop on Heterogeneous Boundary Conditions a general optics problem 
that allows interference was suggested for study. The large-scale interference pattern that 
develops is quite sensitive to small perturbations in the boundary conditions. Hence it 
seemed ideal for testing and observing errors due to grid interface effects introduced by 
domain decomposition methods. Although the problem specification is somewhat arbitrary, 
it is necessary to be specific in order to compare results because it is expected that several 
researchers will explore this problem. 

2 Error Measurements 

The interference that we wish to test lends itself to error analysis using both a visual as well 
as more standard numerical acceptance criteria. The more standard numerical criteria are 

• propagation error 

• L\ error 

• L 2 error 

• Loo error 

By propagation error we mean the phase difference between the computed location of the 
wave front and the exact location of the same wave front. 

Visually this problem gives rise to an interference pattern that can be compared for 
sharpness as well as location. We are interested in seeing these differences across the artificial 
interfaces introduced by the decomposition. 
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3 Motivation 


Dr. Bill Gropp suggested we examine a problem that had visual meaning in its errors, to 
allow studying the types of errors introduced between refined and coarse meshes at a more 
intuitive level. He suggested an optics interference problem. 


4 Optical Interference 

Fermats’ principle of least time is recast in Feynman et. al. [1] briefly from a quantum- 
dynamical perspective. By considering “rays of light” as photons, the ray path can be 
considered a sum of the individual paths of the photons. The ray path is then defined by 
the probability of each photon taking different paths. 

Upon encountering a barrier that contains a wide slit, a wave will continue through 
the slit almost undisturbed and geometric optics is a good model for ray behavior. But 
when the slit is sufficiently reduced in size, the choices for photon paths are truncated, and 
the probability distribution is altered, affecting the geometry of the wave front as it passes 
through the slit. 

This behavior is more easily understood by considering the simpler, less accurate, but 
more intuitive, Huygens principle which states that “all points on a wavefront can be con- 
sidered as point sources for the production of spherical secondary wavelets. After a time the 
new position of the wavefront will be the surface of tangency to these secondary wavelets”, 
as described in Halliday and Resnick [2]. Considering this simplified wave theory of light, a 
barrier in a wave path with a slit on the order of the wave length will cause diffraction. That 
is, the end points allowed to pass through the slit will no longer have symmetrical wavelets 
on either side, and they will bend at the ends. 

Placing two slits aside each other will replicate the limiting behavior twice, and the 
resulting wave patterns will interact, causing interference. The interference pattern will be 
visible, and dependent upon the original wave frequency. This interference pattern is quite 
sensitive to phase errors, so that the choice of grid sizes influences the solution behavior. 
This is the interesting aspect of this problem. 


5 The Interference Pattern 


Consider the wave equation 


d 2 u 

Yt 2 


= Au. 


Here A is the Laplacian. In the positive quadrant place two slits along the x-axis at locations 
X\ and X 2 , with 

d = X 2 — X\. 


Here x* and x 2 are the mid-points of the two slits. Then for any point in the quadrant, 
P = (x p ,y p ), let di be the distance from xj to P, and d 2 be the distance from x 2 to P. 
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Waves will arrive at P out of phase due to the difference in the path lengths d\ and ^ 2 . The 
maximum interference will occur when 

\d 2 — dx\ = m\ 


where A is the original incident plane wavelength and m is a nonnegative integer. The 
minimum, of course, occurs at the half distances (m -f |). 

The size of our slits can now be specified. They should be at most of width A. A slit 
width less than or equal to the incident wavelength is sufficiently small to diffract the wave 
on a visible scale. This problem is interesting because generation of the large scale pattern 
depends on how accurately the small scale dynamics has been captured about the slits. 

* 

6 Geometry 

Here we assume that the problem has been normalized so that 

xi = 0 , 

x r = 1, 

Vb = 0, 

Vt = 1, 

We arbitrarily choose X! = A - 5A, x 2 = | + 5A. This geometry is illustrated on the 
following page in Fig. 1. 

T Initial Conditions 

The initial region, including all boundaries except the slits, should be quiescent (u = 0). 
Prescribe a plane wave described by the following function: 

sin~-(y — vt) (7-1) 

A 

where A is the wavelength and v is the phase velocity. This impacts our domain along the 
slits [A — 5.5A, \ — 4.5A] and [| + 4.5A, | + 5.5A] of the x-axis for all time t > 0. Here, the 
equation reduces to 

O IT 

sin—(—vt). (7-2) 

A 

The phase velocity is known, allowing radiative boundary conditions to be defined. This 
is done in the section 9. 




8 Exact Solution 

The exact interference pattern is the superposition of the two waves. For any point P(x } y) 
in the domain, the travel time to P will be different from the two slits. Let the time from 
slit Xi be Then the value at P will be the sum of the P f where 

Pi{x,y) = sin^-(-v(t - t { )) 

for (l -if) > 0 and zero otherwise. We note that the distance from slit x { to the point P(x, y) 

di = ((x — x ^ 2 + y 2 )* . 


9 Boundary Conditions 

All boundaries, except the slits which are prescribed with the incoming plane wave, evolve 
with the solution. Any workable boundary conditions can be applied on these boundaries, 
with the goal that these boundary conditions should influence the interference pattern as 
little as possible. 

For our prescribed incident plane wave we have the exact solution. However, prescrip- 
tion of these exact boundary values on our numerical approximation of the solution could 
cause the numerical approximation to degrade. Hence we recommend radiative boundary 
conditions because we know the phase velocity exactly and hope that the numerical phase 
velocity is quite close to the correct one. Then an open boundary condition can be used that 
advects the interference pattern out of the domain by advancing the wave equation, 


du du 
dt ^ V dn 


0 


where n is the direction normal to the boundary. 


10 Domain Decomposition Method 

Our purpose in examining this test problem is to measure directly the effect of mesh refine- 
ment and the resulting mesh interfaces on a known wave that is sensitive to phase errors, 
while concurrently being able to visually display a meaningful picture of the effects of the 
refinement-induced error on the solution. 

The Coarse mesh must be able to adequately represent the interference pattern for the 
visual comparisons. It should be no larger than Ax = Ay = and may need to be smaller. 
The workshop suggested We suggest that £, £ and £ all be tested. However, these 
values are only apriori suggestions as we have not yet worked this problem. 

The refined mesh must be able to adequately capture the diffraction behavior, so that 
the plane wave front “bends” as it passes through the slit. Given that the coarse and refined 
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meshes are sufficiently accurate, the phase errors introduced during the problem solution will 
be a function of the sound speed on the two grids plus the coarse/refined grid interaction 
errors. 

Due to the geometry of the problem, non-adaptive local uniform mesh refinement is 
adequate for the domain decomposition. 

11 Discretization 

We suggest a second order in space and time Leap Frog/Hopscotch discretization method 
should be used. We anticipate beginning with A = 0.1, and A = 0.05. 

12 Analysis 

We expect that calculations for boundary errors and the ratio of mesh refinements will be 
analyzed, These can be done both analytically and by comparison to an everywhere fine 
mesh. 
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